#include "Numeric.h"

using namespace std;
//***** Scott Butler's code Start
//***member functions of class matrix
Matrix::Matrix() {}
Matrix::Matrix (int rows, int cols) {
    this->rows=rows;
    this->cols=cols;
    this->data.resize(rows*cols);
    for(int r=0; r<rows; r++)
      for(int c=0; c<cols; c++)
	this->data[r*cols+c]=0;
  }

double Matrix::operator()(int i, int j) const {
    return data[i*cols+j];
  }

double & Matrix::operator()(int i, int j) {
    return data[i*cols+j];
  }

Matrix Matrix::operator+(const Matrix &B) {	
  Matrix A(rows,cols);//contruct calling matrix
  for(int r=0; r<rows; r++)//copies calling matrix
   for(int c=0; c<cols; c++)
      A(r,c)=data[r*cols+c];
  if(A.rows!=B.rows || A.cols!=B.cols)
    cout << "BAD\n";
  Matrix C(A.rows,A.cols);
  for(int r=0; r<B.rows; r++)
    for(int c=0; c<B.cols; c++)
      C(r,c)=A(r,c)+B(r,c);
  return C;
}

Matrix Matrix::operator-(const Matrix &B) {
  Matrix A(rows,cols);//contruct calling matrix
  for(int r=0; r<rows; r++)//copies calling matrix
   for(int c=0; c<cols; c++)
      A(r,c)=data[r*cols+c];
  if(A.rows!=B.rows || A.cols!=B.cols)
    cout << "BAD\n";
  Matrix C(A.rows,A.cols);
  for(int r=0; r<A.rows; r++)
    for(int c=0; c<A.cols; c++)
      C(r,c)=A(r,c)-B(r,c);
  return C;
}

Matrix Matrix::operator-(const double &x) {
  Matrix A(rows,cols);//contruct calling matrix
  for(int r=0; r<rows; r++)//copies calling matrix
   for(int c=0; c<cols; c++)
      A(r,c)=data[r*cols+c];
  Matrix C(A.rows,A.cols);
  for(int r=0; r<A.rows; r++)
    for(int c=0; c<A.cols; c++)
      C(r,c)=A(r,c) - x;
  return C;
}

Matrix Matrix::operator*(const Matrix &B) {
  Matrix A(rows,cols);//contruct calling matrix
  for(int r=0; r<rows; r++)//copies calling matrix
   for(int c=0; c<cols; c++)
      A(r,c)=data[r*cols+c];
  if(A.cols!=B.rows)
    cout << "BAD\n";
  Matrix C(A.rows,B.cols);
  for(int r=0; r<A.rows; r++)
    for(int c=0; c<B.cols; c++)
      for(int k=0; k<A.cols; k++)
	C(r,c)+=A(r,k)*B(k,c);
  return C;
}

Matrix Matrix::operator*(const double &a) {//must mutiply Matrix first * scalar second
  Matrix A(rows,cols);//contruct calling matrix
  for(int r=0; r<rows; r++)//copies calling matrix
   for(int c=0; c<cols; c++)
      A(r,c)=data[r*cols+c];
  Matrix C(rows,cols);
  for(int r=0; r<A.rows; r++)
    for(int c=0; c<A.cols; c++)
      C(r,c)=a*A(r,c);
  return C;
}


void Matrix::swap(double&a, double &b) {
  double c=a; a=b; b=c;
}

Matrix Matrix::inv(Matrix A) {
  if(A.cols!=A.rows)
    cout << "BAD\n";
  Matrix B(A.rows,A.cols);
  double p;
  double q;
  int m;
  for(int r=0; r<B.cols;r++) B(r,r)=1;
  for(int c=0; c<A.cols;c++) {    
    m=c; p=A(c,c);
    for(int i=c+1; i<A.rows; i++)
      if(abs(A(i,c)) > abs(p)) {m=i; p=A(i,c);}
    for(int i=0; i<A.cols; i++) {
      swap(A(m,i),A(c,i));
      swap(B(m,i),B(c,i));
    }
    for(int i=0; i<A.cols; i++) {
      A(c,i) /= p; 
      B(c,i) /= p;
    }
    for(int r=0; r<A.rows; r++) 
      if(r!=c) {
	q = A(r,c);
	for(int i=0; i<A.cols; i++) {
	  A(r,i)-=q*A(c,i);
	  B(r,i)-=q*B(c,i);
	}
      }
  }
  return B;
}

Matrix Matrix::identity(int i){
	Matrix A(i,i);
	for(int r=0; r<A.rows; r++)
		for(int c=0; c<A.cols; c++)
			A(r,c)= 0;
	for (int r = 0; r<A.rows; r++)
		A(r,r) = 1;
	return A;
}
Matrix Matrix::operator/(const double &scalar) {//matrix/scalar
	Matrix A(rows,cols);//contruct calling matrix
	  for(int r=0; r<rows; r++)//copies calling matrix
	   for(int c=0; c<cols; c++)
		  A(r,c)=data[r*cols+c];
	  return A * (1/scalar);
}
Matrix Matrix::operator /(const Matrix &B) {//matrix/matrix
	Matrix A(rows,cols);//contruct calling matrix
	  for(int r=0; r<rows; r++)//copies calling matrix
	   for(int c=0; c<cols; c++)
		  A(r,c)=data[r*cols+c];
	return A * (1.0/B);
}
Matrix operator/(const double &scalar, Matrix A) {//gauss jordan
	int n = A.cols;
	int indexes = A.cols;
	if (n !=A.rows) cout << "Matrix not squared" << endl;
	Matrix B(n, n);
	B = B.identity(n);
	B =  B * scalar;
	for (int c = 0; c < indexes; c++) {
		for (int r = c+1; r < n; r++) {
			if (abs(A(r,c)) > abs(A(c,c))){
				for (int swapCol = 0; swapCol < A.cols; swapCol++){
					swap((A(r, swapCol)),A(c, swapCol));
				    swap((B(r, swapCol)),B(c, swapCol));
				}
			}
		}
		double p = A(c,c);
		for (int k = 0; k < indexes; k++) {
			A(c,k) = A(c,k)/p;
			B(c,k) = B(c,k)/p;
		}
		for(int r= 0; r < n; r++) if(r!=c){ 
			p = A(r,c);
			for(int k = 0; k < indexes; k++) {
			A(r,k) -= A(c,k) * p;
			B(r,k) -= B(c,k) * p;
				}
		}
		
	}
	return B;
}
Matrix Matrix::tranpose(const Matrix &A) {
	Matrix newMatrix(cols,rows);
	for (int r = 0; r < rows; r++)
		for (int c = 0; c < cols; c++) 
			newMatrix(c,r) = A(r,c);
	return newMatrix;
}
Markowitz::Markowitz() {
portReturn = 0;
portRisk = 0;
}

//input R as return vector, C as cov matrix of stock returns, and risk free rate as risk free
//loads optimize pf weights into portReturns, port
void Markowitz::compute_markowitz(Matrix returnVector, Matrix covMatrix, double riskFree) {
	double sum=0;
	Matrix invRealCovMatrix(covMatrix.rows,covMatrix.cols), risk(1,1);

	//computation
	invRealCovMatrix = (1/covMatrix) * (returnVector-riskFree);
	for (int i=0; i<invRealCovMatrix.rows;i++) sum+=invRealCovMatrix(i,0);
	invRealCovMatrix = invRealCovMatrix/sum;

	//load into class private variables that are portflio weights 
	portWeight.resize(invRealCovMatrix.rows);
	for (int i=0; i<invRealCovMatrix.rows;i++) {
		portWeight[i] = invRealCovMatrix(i,0);
		portReturn += (invRealCovMatrix(i,0) * returnVector(i,0));
	}
	
	risk = (invRealCovMatrix.tranpose(invRealCovMatrix) * (covMatrix * invRealCovMatrix));
	portRisk = sqrt(risk(0,0));
}

//no input.  returns matrix of opitmized portfolio weights after compute_markowitz is run
Matrix Markowitz::optimized_portfolio_wts(){
	Matrix portWeightMatrix(portWeight.size(), 1);

	for (int i=0; i<portWeightMatrix.rows; i++) portWeightMatrix(i,0) = portWeight[i];
	return portWeightMatrix;
}

double Markowitz::optimized_portfolio_rtn() {
	return portReturn;
}

double Markowitz::optimized_portfolio_risk() {
	return portRisk;
}
double FunctionScott::D(double x, double h=1e-6) {
    return (f(x+h)-f(x-h))/(2.0*h);
}


//Norm Function, must input double vector, matrix, or a double/int value
//vector.
double norm(const vector<double> &x) {
	double sum = 0, p=1;
	for (int i = 0; i < x.size(); i++)
		sum += (abs(pow(x[i],p))*1.0/p);
	return sum;
}
//normal for matrix, input matrix returns norm of matrix
double norm(const Matrix &A) {
	vector<double> elementHolder;
	double sum = 0, p=1;
	if (A.rows==1) {
		for (int i=0; i<A.cols;i++) sum += abs(A(0,i));
		return sum;
	}
	if (p==1) {
		for (int c = 0; c < A.cols; c++){
			for (int r = 0; r < A.rows; r++)
				sum += abs(A(r,c));
			elementHolder.push_back(sum);
			sum = 0;
		}
		return *max_element(elementHolder.begin(),elementHolder.end());
	}
	else cout << "not implemented";
}
//norm for double, returns abs value
double norm(const double &x) {
	return abs(x);
}

//symmetric function.  input matrix.  returns true or false is matrix is symmetric
bool is_almost_symmetric(const Matrix &A) { 
	double ap=1e-6, rp=1e-4, delta;
	if(A.rows != A.cols) return false;
	for (int r = 0; r < A.rows; r++) {
		for (int c=0; c < r; c++) {
			delta = abs(A(r,c) - A(c,r));
			if(delta > ap && delta > (max(abs(A(r,c)), abs(A(c,r)))*rp)) 
				return false;
		}
	}
	return true;
}

//input matrix .  returns true or falase if matrix is almost zero
bool is_almost_zero(const Matrix &A) {
	double ap=1e-6, rp=1e-4, delta;
	for (int r = 0; r < A.rows; r++) {
		for (int c=0; c < r; c++) {
			delta = abs(A(r,c) - A(c,r));
			if(delta > ap && delta > (max(abs(A(r,c)), abs(A(c,r)))*rp)) 
				return false;
		}
	}
	return true;
}

// condition number imput matrix, returns condition number
double Condition_Number::matrix_cn(const Matrix &A) {
	double norm1, norm2;
	norm1 = norm(A);
	norm2 = norm(1.0/A);
	return norm1*norm2;
}
double Condition_Number::func_cn(double x) {
	return  D(x)*x/f(x);
}
double Condition_Number::f(double x){//just to test 
	return x*x*x;
}

//input matrix and returns a matrix after choleksy is performed. 
//if not symmetric or pos defiinite function aborts
Matrix cholesky(const Matrix &A) {
	double p;
	Matrix L(A.rows,A.cols);
	L = A;
	if(!is_almost_symmetric(A)) {
		throw string("Not Symmetric");
	}
	for (int k = 0; k < L.cols; k++) {
		if (L(k,k) <= 0) {
			throw string("Error:  Not Positive Definite");
		}
		p = sqrt(L(k,k));
		L(k,k) = p;
		for (int i =k+1; i<L.rows; i++) {
			L(i,k) /=p;
		}
		for (int j = k+1; j<L.rows; j++){
			p = L(j,k);
			for (int l = k+1; l<L.rows; l++) 
				L(l,j) -= p*L(l,k);
		}
	}
	for (int i=0; i < L.rows; i++) {
		for (int j=i+1; j <L.cols; j++)
			L(i,j) = 0;
	}
	return L;		
}

bool is_positive_definite(const Matrix &A) {
	if(!is_almost_symmetric(A)) return false;
	try {
		cholesky(A);
		return true;
	}
	catch(...) {
		return false;
	}
	return true;
}
ostream &operator <<(ostream &out, const Matrix& A) {
  out << "[";
  for(int r=0; r<A.rows; r++) {
    if(r>0) out << ",";
    out << "[";
    for(int c=0; c<A.cols; c++) {
      if(c>0) out << ",";
      out << A(r,c);    
    }
    out << "]";
  }
  out << "]";    
  return out;
}
//***** Scott Butler's code End